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A  new  method  for  reconstruction  of  the  shape  of  the  left  ventricle  from  the 
biplane  angiocardiograms  is  proposed.  This  approach  utilizes  a  pair  of  orthogonal  X- 
ray  projection  images.  The  shape  of  the  ventricle  is  reconstructed  by  dividing  these 
projection  images  into  parallel  slices  and  then  being  processed  slice  by  slice 
stepwise.  Each  corresponding  pair  of  slices  form  two  one-dimensional  projection 
profiles  which  are  used  to  reconstruct  a  cross  section  of  the  ventricle. 

Without  using  predefined  models,  we  proposed  a  new  method  to  reconstruct  the 
cross  section  under  the  assumption  that  the  cross  section  is  regular  along  the 
directions  of  projections  and  a  monotonically  nondecreasing  or  nonincreasing  equal- 
divisor  curve  is  available.  A  cross  section  is  regular  if  it  contains  only  one  closed 
interval  within  the  cross  section  along  each  ray  of  projection.  The  equal-divisor 
curve  is  defined  within  the  cross  section  which  divides  each  closed  interval  into  two 
equal  halves.  Instead  of  solving  the  binary  matrix  directly,  we  need  only  to  obtain 
the  equal-divisor  curve,  then  the  cross  section  is  uniquely  determined-  We  have 
proved  that  the  cross  section  can  be  optimally  reconstructed  by'  minimizing  an  error 
index  which  is  generated  by  the  sum  of  the  absolute  difference  between  the  original 
and  estimated  projection  profiles.  The  algorithm  has  been  tested  on  synthetic  and  X- 
ray  pictures  with  better  than  95%  conformity  for  regular  cross  sections.  The 
experiments  show  better  results  than  any  existing  direct  methods  since  it  is  based 
on  a  more  flexible  geometrical  assumption.  It  also  prevents  the  sensitive  selection  of 
mask  model  which  dominates  the  results  of  model-based  methods  Thus,  the 
algorithm  is  reliable,  and  more  generally  fits  real  shapes  of  the  ventricle 


CHAPTER  1 


INTRODUCTION 

In  patients  with  cardiovascular  diseases,  left  ventricular  function  is  one  of  the 
major  parameters  to  be  considered  in  making  therapeutic  choices  and  determining 
prognosis.  Left  ventricular  function  is  primarily  evaluated  by  measurements  of 
pressure  and  volume  during  a  cardiac  cycle.  Although  several  indices  based  upon  the 
volume  measurements  have  been  developed,  accurate  measurements  of  volume  are 
not  yet  possible.  Almost  all  methods  of  measurement  are  based  upon  an  ellipse 
model.  Even  the  biplane  integration  method  based  on  Simpsons  rule,  the  most 
accurate  method  of  all.  assumes  an  ellipsoid  shape  constructed  around  two  axes  for 
each  cross  section.  Therefore,  any  method  that  can  improve  the  accuracy  of  the 
measurement  of  cross  sectional  area  is  likely  to  improve  accuracy  of  the  method  to 
measure  the  ventricular  volumes.  The  present  study  was  undertaken  to  investigate  if 
the  shape  of  each  cross  section  can  be  predicted  with  better  accuracy  knowing  the 
density  of  the  dye  (radiopaque  dye  injected  to  visualize  the  left  ventricle)  in  two 
planes. 

In  this  report  we  present  a  new  method  for  reconstruction  of  the  shape  of  the 
left  ventricle  from  the  biplane  angiocardiograms.  This  approach  utilizes  a  pair  of 
orthogonal  X-ray  projection  images  for  each  ventricular  phase  in  a  cardiac  cycle 
The  projection  image  is  a  two-dimensional  picture  which  is  obtained  by  integrating 
the  absorption  function  of  the  X-ray  photons  along  the  directions  parallel  to  the 
optical  axis  of  the  camera  sets.  Obviously,  the  absorption  function  affected  by  the 
injected  dye  will  influence  the  reconstructed  results  For  simplicity,  we  define  the 
absorption  coefficient  of  the  dye  to  be  unity  to  form  a  binary  reconstruction 
problem  just  like  ail  the  other  developed  methods  The  shape  of  the  ventricle  is 
reconstructed  by  dividing  these  projections  into  slices  and  then  being  processed 


slice  by  slice  stepwise  Therefore,  we  reduce  our  three-dimensional  problem  into  a 
two-dimensional  problem  by  considering  the  object  as  a  stack  of  parallel  cross 
sections.  Each  cross  section  consists  of  one  connected  region  and  is  reconstructed 
from  its  two  one-dimensional  projection  profiles.  Therefore,  our  problem  can  be 
stated  as  the  reconstruction  of  a  binary  matrix  from  its  row  and  column  sums. 

In  Chapter  2  we  describe  our  reconstruction  problem  and  review  some  currently 
developed  methods.  In  Chapter  3  some  mathematical  characteristics  are  explicated. 
Our  algorithm  and  its  implementation  are  discussed  in  Chapter  4.  Some  synthetic 
examples  are  also  demonstrated.  In  Chapter  5,  we  describe  an  experiment  by  using 
two  projection  images  from  a  bag  of  dye.  The  discussion  of  the  experimental 
results  and  some  suggestions  for  further  improvements  are  given  in  Chapter  6 


CHAPTER  2 


PROBLEM  DESCRIPTION  AND  LITERATURE  REVIEW 

The  binary  reconstruction  problem  may  be  stated  as  follows  in  terms  of  entries  of 
a  matrix  X  with  m  x  n  cells.  Each  cell  is  identified  by  the  pair  of  numbers  (i,j>  and  is 
represented  by  x ...  x  may  be  1  if  the  cell  falls  into  the  region  of  object  or  0 
otherwise.  The  problem  is  to  estimate  x  given  the  marginals 

n 

l  x. .«P. ,  i «1 ,  ...  ,m 

j=i 

m 

l  xm“V  J-1*  •  •  *n 

i«1  J  J 

where  P  i=1,  .  .  ,m  are  elements  of  vector  P  containing  the  row  sum  of  X,  and 
Q  j=1,..,n  are  elements  of  vector  Q  containing  the  column  sum  of  X.  It  is  also  clear 
from  the  essence  of  the  matrix  that 

m  n 

I  p,-Z  Q; 

i*  1  j=1  J 

This  implies  that  the  estimation  of  m  x  n  variables  x  e{0, 1}  needs  to  be 

'j 

determined  from  m  +  n  -  1  independent  equations.  Therefore,  our  problem  is 
obviously  underdetermined.  Some  previous  research  results  (2,5,7)  also  show  that 
the  reconstruction  of  an  arbitrary  binary  matrix  from  two  projections  is  in  general 
impossible.  According  to  the  theory  by  Ryser  (1,2),  it  is  possible  to  find  zero,  one, 
or  more  than  one  binary  matrices  satisfying  a  given  pair  of  row  and  column  sums 
(see  Fig.2-1).  To  reconstruct  the  binary  matrix,  some  certain  assumptions  always 
need  to  be  made.  In  the  following,  we  are  going  to  review  some  of  the  research 
works  done  in  this  area  and  their  assumptions 
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The  matrix  reconstruction  algorithms  can  roughly  be  divided  into  two  categories 
The  algorithms  in  the  first  category  apply  some  direct  methods  to  solve  the 
matrices  under  some  assumptions,  while  those  in  the  other  category  usually 
reconstruct  the  matrix  by  using  a  priori  knowledge,  in  other  words,  some 
predefined  binary  mask  models.  We  will  discuss  some  methods  in  both  categories  as 
follows. 

The  commonly  used  method  in  biplane  angiocardiograms  utilizes  the  assumption  that 
the  shape  of  the  cross  section  is  close  to  an  ellipse  (3).  So  the  distance  between 
two  end  points  of  each  projection  profile  is  taken  as  an  axis  of  the  ellipse 
However,  the  shape  of  the  cross  section  is  usually  not  an  ellipse,  and  the 
reconstruction  result  of  a  certain  cross  section  also  depends  on  the  orientation  of 
the  cross  section.  In  1983,  Eiho  et  al.  (4)  proposed  a  method  to  solve  the 
orientation  dependence  problem  by  using  three  projection  profiles  with  the  same 
ellipse  assumption.  The  reconstructed  shape  is  much  improved  with  the  expense  of 
one  more  projection,  but  still  far  from  a  shape  of  a  real  ventricular  cross  section. 
Chang  (5)  has  developed  an  algorithm  to  recover  any  binary  matrix  from  two 
projection  profiles.  The  algorithm  can  obtain  the  exact  reconstruction  when  there 
exists  only  one  unambiguous  solution  from  two  projection  profiles  However,  most 
of  the  shapes  of  cross  sections  are  ambiguous  based  on  his  definition.  Only  one  in 
a  large  amount  of  possible  solutions  found  by  his  algorithm  is  the  real  solution  So 
the  results  are  usually  not  correct  Chang  and  Chow  proposed  another  algorithm  to 
specifically  reconstruct  the  shape  of  a  ventricle  (6).  They  assumed  the  cross  section 
of  the  ventricle  should  be  in  a  connected  region  which  is  convex  symmetric  So  the 
ambiguity  is  largely  reduced,  only  a  few  possible  solutions  remain  after 
reconstruction.  However,  a  cross  section  of  the  left  ventricle  is  in  general  not 
convex  symmetric  (7,9).  Therefore,  we  proposed  our  new  algorithm  without  the 
convex  symmetry  assumption  to  fit  more  generally  to  the  left  ventricle 

Because  the  ambiguity  was  hardly  solved  by  most  of  the  direct  methods, 
researchers  also  found  another  approach  that  uses  a  priori  information.  At  first, 
Onnasch  and  Heintzen  (7)  and  Onnash  (8)  developed  an  algorithm  for  binary 
reconstructiobn  of  the  left  or  right  ventricle  from  two  projections  that  utilizes  a 


mode!  database  obtained  from  the  cast  study  of  ventricles.  This  method  implies  the 
estimation,  for  every  matrix  element,  of  its  probability  to  belong  to  the  ventricle. 
Recently,  Slump  and  Gervranchs  proposed  a  method  also  incorporated  with  a  priori 
knowledge  to  reduce  the  ambiguity  of  the  problem  (9).  A  minimum  cost  capacitated 
network  flow  algorithm  is  adopted,  which  yields  the  optimal  solution  with  respect  to 
the  selected  binary  matrix  models.  If  the  correct  models  are  available,  these  model- 
based  reconstruction  methods  usually  give  more  promising  results  than  the  direct 
methods.  However,  a  correct  model  is  physically  unavailable,  in  fact,  it  is  the  goal  to 
be  recovered  by  these  algorithms  Our  new  algorithm  reconstructs  the  binary  matrix 
with  the  assumption  that  the  equal-divisor  curve  of  a  cross  section  always  exists 
No  binary  matrix  model  is  required,  though  the  reconstruction  results  are  pretty 
reliable.  We  will  discuss  it  in  the  following  sections. 


CHAPTER  3 


MATHEMATICAL  CONSIDERATIONS 


Suppose  that  V  is  a  connected  region  in  R  (d  >=  21  and  fix)  =  1  v<x),  the  indicator 
function  of  V.  Define 

p(x1 - Vl1^  f  (W  '  '  ’Xd)dXd 


and 


q  (x  ,  .  .  ,x  )=/  f  (x.,x.,  .  .  .  ,x  )dx 
2  a  l  2  a  I 


(3  1) 


(3.2) 


p  and  q  are  known  as  projections  along  x  -axis  and  x  -axis  respectively.  In 

\  d 

practice,  our  observations  are  given  as  follows 

R1  (x  i ,  .  .  ,xd_1)  «=p  (x, ,  .  .  ,xd_1)  +  e(x1,  .  .  •xd.1)+X1 


W 


!q  (x. 


x  .)+ti  (x  , 

Q  2 


Xd)+X2 


(3-3) 


(3.M 


where  e  and  ti  are  known  as  noises  centered  at  expectations  and  Ai  .  X2  are 
known  as  the  expectations  of  noises.  In  continuous  case,  e  and  ti  can  be  assumed 
to  be  stochastic  processes  and  in  descrete  case,  they  are  assumed  to  be  two 
families  of  random  variables  In  practice.  Ri  and  can  be  computed  from  two 
photos  which  are  taken  along  two  orthogonal  directions.  The  problem  we  are  faced 
is  how  to  reconstruct  the  shape  of  the  subject  V  by  using  these  two  profiles.  R 
and  R  The  mathematical  work  concerning  this  problem  includes 

(i)  estimation  of  the  projections  p  and  q 


g 


(ii)  reconstruction  of  V,  V,  by  using  p  and  q.  < 

i 

l 

Obviously,  the  first  task  is  statistical  and  the  second  is  geometrical.  In  this  section,  j 

we  shall  discuss  how  to  solve  these  two  problems.  < 

t 

t 

t 

At  first,  we  shall  consider  the  noise  free  case.  In  this  case,  there  is  no  need  to  J 

estimate  the  projections  and  it  is  assumed  that  p  and  q  are  known,  A1  and  X2  are  \ 

zeros.  In  our  reconstruction  approach,  the  two  projection  profiles  are  not  equally 
treated.  According  to  their  usage,  let  us  name  p  as  the  constructor  (constructive 
profile)  and  q  as  the  ruler  (control  profile),  respectively.  j 

( 

In  the  sequel,  we  always  assume  that  V  are  regular,  namely,  each  straight  line 
along  x  -axis  or  x  -axis  intersects  V  in  a  closed  interval  on  the  line,  if  they  are  not 

1  d 

disjoint  Define  the  equal-divisor  surface  (ED)  as  follows: 


D  =  {(x  ,  .  .  .  ,x  ),  there  is  an  X  such  that  (x . x  )  eV} 

1  d- 1  d  id 


For  (x  ,  .  .  .  ,x  )  e  D, 

1  d-1 

denote  the  two 

intersection 

points 

of  the  line 

lxrx, . X,.,=V,'  and 

(X  .X . x'l,  i.e., 

1  2  d 

the  surface  of 

V  by  (x 

TX2' 

°v 

.  .  ,x  )  and 
d 

x°  =  sup  {x  ,  (x  ,x  ,  . 
d  d  12 

•  ,x  )  eV} 

d 

and 

x1  =  inf  {x  ,  (x  ,x  ,  . 
d  d  12 

.  .  ,x)  eV} 

d 

Then  the  surface  {(x  ,  .  .  .  ,x 

1  d 

,(x°  +  x  V2  )  ,  (x 

—  Id  d 

. x  ) 

1  d-i 

1  e  D} 

is  called  the 

ED  of  V  (or  more  definitely  the  ED  of  V  along  x  -axis).  For  convenience,  we  also 

d 

call  the  function  g(x  ,  .  .  .  ,x  )  =  l/2(x  °+  x  ')  the  ED  of  V. 

1  d- 1  d  d 


The  following  propositions  are  evident. 


Proposition  1.  If  the  ED  and  the  constructor  of  V  are  known,  then  V  can  be 
reconstructed  uniquely. 


10 


Proposition  2.  If  the  constructor  is  given,  then  V  and  its  ruler  are  continuous  with 
respect  to  ED  in  the  following  metrics: 
d  (V,V)  -/  |  lv-lv,  |dx 

(3-5) 


d(q.q')-/  |  q  (x  ^  •  •  .  xd>  ~q  *  (x2>  ,xrf)  |dx2>  ,dxd 


(3-6) 


d(g.g')-/  |g(xr  .  . xd_ ^ ~s '  (*,.  ■  .x  j)  Idx^  .  ,dxdJ 


where  g  =  1/2(x  °+x  S  and  g'  is  similarly  defined  for  another  subject 

d  d 


(3-7) 


The  most  difficult  thing  in  reconstruction  is  that  there  is  no  unique  solution  We 
have  the  following  example. 

Example  1  Suppose  and  are  two  identical  ellipses  except  the  long  axis  of 

Vi  coincides  with  the  line  y=x  and  the  long  axis  of  V2  coincides  with  y=-x.  Both 

of  their  centers  are  at  origin.  Then  V  and  V  have  the  same  constructors  and 

1  2 

rulers 

To  make  the  solution  unique,  we  have  to  make  further  assumptions  which  can  be 
justified  by  practical  knowledge  For  instance,  the  doctor  can  say  the  heart  of  a 
patient  is  tilted  towards  left  by  his  experience 

In  the  following  we  will  consider  the  case  d  =  2  Suppose  p(x|  is  a  continuous  and 
unimodal  function  defined  on  the  interval  D  =  [a,b]  with  p(a)  >  0.  p(b)  >  0.  and  p(x) 
>  0.  x  e(a.b).  Let  f  be  a  function  defined  on  [a,b]. 

Define 

V  (f)-{  (x,y)  :xeD,f  (x)  -l/2p(x)  <y<f  (x)  +  l/2p  (x) } 

and  define 

q  ^  (y)“hax{x)-x2;  (x]  ,y)  eV  (f)and  (x^y)  r  (e)  V(f)  } , 

(3-8) 

if  V(f)  is  regular  and  if  there  is  at  least  one  point  x  such  That  (x.y)  e  V(f),  and 
q<fl(y)  =  0  otherwise 


1 1 

Suppose  F  is  a  class  of  functions  defined  on  [a,b]  satisfying: 

(i)  each  f  in  F  is  nondecreasing  and  continuous, 

(ii)  for  each  f  eF,  V(f)  is  regular, 

(iii)  for  any  pair  f  ,  f  t  F,  the  set  {x  e  (a,b),  f  (x)  P  f  (x)}  consists  of  at  most 
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two  pieces  of  open  intervals, 
v  =  {  V(f)  :  f  e  F} 

(3.9) 

and 

R  -  {q(f)  : f  e  F} 

(3. 1 0) 

We  have  the  following  basic  theorem. 

Theorem  3.1  Under  the  assumptions  given  above,  we  have 
F  V  R 

(3.1  1) 

Where  denotes  one-to-one  correspondence.  The  proof  of  Theorem  3.1  is  given 
in  a  companian  paper  by  some  of  the  present  authors. 

Theorem  3.2  (Basic  Reconstruction  Theorem).  Suppose  V  is  a  subject  to  be 
reconstructed,  and  suppose  the  constructor  of  V  is  P(x).  Let  F  be  a  function  class 
satisfying  the  conditions  given  in  (3.9).  Suppose  the  ED  of  V  belongs  to  F  (of 
course,  this  implies  V  is  regular).  Then  V  is  uniquely  determined  by  its  ruler  among 
V  with  respect  to  F. 

2 

Example:  Suppose  V  is  convex  region  in  R  with  ED  being  a  straight  line  up- 
slooed  P(x),  x  e  [0,a]  and  q(y),  y  €  [0,b]  are  its  constructor  and  ruler 
respectively.  Take 

F  =  {ax  +  3-  0  <  a  <  b/a.  0  <  (3  <  b} 

Use  P  and  F,  we  get  the  set  V,  then  by  Theorem  3.2  V  is  the  only  element  in  V 


which  has  q(y)  as  its  ruler.  This  example  explains  why  our  procedure  given  in 

section  4-1  works.  Now  let  us  extend  the  result  from  d  =  2  to  the  general  case. 

Let  V  be  a  convex  region  in  Rd,  d  >  2.  As  earlier  defined,  let  p  and  q  be  its  two 

projections  and  let  g  be  its  ED  For  any  fixed  x  =  x  0 . x  =  x°  ,  define 

9  22  d-1  d- 1 

*  •  0  0 
v  “C(x,»xd)  :  (x1  ,x  ,  ...  » Xrf_  1 . xd)  E  V] 

*  0  0 

P  (x^-p  (x,  ,x2,  .  .xd1) 

q*(xd)~q(x°,  .  .x°_1,xd) 

and 

*  0  0 

g  (x^-g  (x]  ,x2,  •  ■  ,xd  l) 

ft  ft 

provided  that  V  is  not  empty.  Then  V  can  be  regarded  as  the  cross  section  of  V, 

0  0  *  #  # 

cut  out  by  the  hyper  plane  {x  =x  . x  =  x  }.  Also,  p  ,  q  and  g  are  the 

22  d-1  d-1  '  ' 

ft 

constructor,  ruler  and  ED  of  V  respectively.  It  is  easy  to  see  that  for  each 
o  o  * 

x2=x2 . Xd-i=X  d  i'  9  'S  continuous  and  nondecreasing  if  and  only  if  g  is 

continuous  and  nondecreasing  in  x  for  each  fixed  x . x  .  Hence  by  theorem 

1  2  d- 1 

3.2  we  have 

Theorem  3.3  Suppose  that  V  is  a  regular  subject  in  Rd  to  be  reconstructed  with 

nondamaged  profiles  P(xi . Xd- ^  anci  . Xd^‘  Let  F  be  a  family  of 

functions  satisfying  the  conditions  (3.9)  for  each  given  multiple  (x  ,  .  .  .  ,x  ).  Then 
the  true  subject  V  is  the  only  one  which  has  q  as  its  ruler,  among  the  family  of 
subjects  constructed  by  p  and  each  element  in  F. 

d  2 

Here  the  regularity  of  region  in  R  is  similarly  defined  as  that  in  R  .  Namely,  a 

region  V  in  Rd  is  called  regular,  if  each  line  {x  =  x  0 . x  =x°  }  or 

9  1  1  d- 1  d- 1 

r  0  0, 

lx2=x2 . x^=  x^  j  can  only  intersect  V  in  a  closed  interval  (or  a  single  point 

or  empty). 

Remark  1.  Since  the  regions  under  usual  consideration  are  assumed  to  be  convex, 
thus  the  condition  of  regularity  is  satisfied. 
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By  Theorems  3.2  and  3.3,  the  task  of  reconstruction  turns  out  to  estimate  the  ED 
of  the  region  by  using  its  two  projections.  Thus  our  problem  becomes  to  fine  a 
function  among  a  given  class  of  functions  such  that  the  derived  projection  along 
x  -axis  is  closest  to  the  observed  one,  the  ruler. 

i 

By  Weierstrassian  approximation  theorem,  for  any  continuous  function,  we  can  find 
a  polynomial  uniformly  approximating  it  Thus  we  can  restrict  our  attention  to  finding 
a  polynomial  as  the  estimate  of  the  ED  of  V.  The  detail  description  of  the 
reconstruction  algorithm  will  be  given  in  the  following  section. 

Next,  let  us  consider  the  case  where  noise  arises.  Because  the  projections  p  and 
q  are  damaged  by  the  noise,  we  can  not  use  and  to  do  reconstruction  and 
we  have  to  smooth  them  to  eliminate  the  affect  of  noise.  There  are  many  ways  to 
do  so.  We  only  state  several  simple  methods  here. 

(i)  Kernel  smoothing.  Suppose  K  is  a  probability  density.  Choose  a  positive  number 
h.  then  let 

P  (x)  +Xi*  -  /  R,(y)K( — ~)  dy 

hd_1  Rd-’  H 

q(x)+X,*=  -  /  R_  (y ' )  K  ( - — }dy' 

hd_1  „<*-!  1  h 

h  r 

(ii)  Histogram  smoothing.  Split  D  into  small  regions  as  D  ,  D . D  .  Then 

1  2  m 

define 

p(x)+X1*  —  -  j  /  R  j  (y ' )  dy ,  i  f  x  e  D. 

Similarly  smooth  q  as  q. 

(iii)  Median  smoothing.  Choose  6  >  0,  for  each  x,  take  the  median  of  the 

A 

observations  of  p  observed  in  the  interval  (x~6,x+6)  as  the  estimate  p  of  R  . 

a 

Similarly  define  q.  This  approach  is  available  only  for  discrete  case. 

(iv)  Polynomial  fitting.  Find  a  polynomial  p  which  minimizes 
/  [R,  (x)  -  (p(x)l^  ^  (x)  +X)  ]2dx 

A 

Similarly  define  q. 


Under  certain  conditions,  we  can  prove  all  the  above  approaches  get  consistent 
estimates  of  the  true  profiles  The  proof  will  be  found  on  the  other  technical  report 
which  will  be  filled  later  soon. 


CHAPTER  4 


SHAPE  RECONSTRUCTION 


For  a  given  cross  section  as  in  Fig.  2,  one  can  have  the  two  profiles  from  the 
column  and  row  projections,  where 


m 

I 

i=1 

X.  . 

IJ 

•  P  j .  j  -  1  , 

.n  ; 

m 

I 

i=i 

X.  . 

U 

i 

O’ 

a 

.  .  ,m  ; 

(4.1) 


P  and  Q  denote  the  two  projection  profiles  and  x  's  e  {0.1}  are  the  cells  in 

i  '  'j 

the  given  matrix.  We  assume  that  a  cross  section  of  a  left  ventricle  is  a  connected 
region  without  holes  in  the  m  x  n  binary  matrix  X  Therefore,  we  can  get  the  mass 
center  of  the  cross  section  which  can  always  be  unambiguously  determined  from 
any  given  sets  of  profiles. 


i 

c 


m 

I  '  *Q; 

i*  1 

-  and 

m 

i «. 

i=  i 


l  i 


1*1 


p. 

j 


(4.2) 


In  general,  we  can  find  an  ED  curve,  along  which  one  of  the  projection  profiles 
can  be  divided  into  two  equivalent  halves  For  instance,  in  Figure  4-1,  we  can  see 
the  real  ED  curve  divides  P  into  two  halves  at  j  =  j,  where  one  half  is  above  the 

j 

ED  curve  and  the  other  half  beneath  it  Consequently,  if  we  repeat  this  operation 
for  j  =  0,  .  ,n.  along  the  real  ED  curve,  the  cross  section  can  be  correctly 
reconstructed  from  P  In  other  words,  the  proposition,  that  a  cross  section  can 

i 
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always  be  uniquely  determined  by  one  of  its  projection  profiles  and  its 
corresponding  ED  curve,  is  always  valid  for  any  regular  cross  sections  Therefore, 
our  problem  is  simply  to  estimate  the  real  ED  „jrve  with  respect  to  P  or,  similarly 
the  real  ED  curve  with  respect  to  Q,  which  can  divide  Q,  i  =  1,  ,m  into  two 

l  l 

equivalent  halves  rowwise  Let  us  use  the  following  equation  to  approximate  the  ED 
curve  in  Figure  4-1. 

i  ( j ' )  »  ic+a0+a1  (j'~jc)+a2(j  '~jc)2+a3  (j'"jc)  3 

(4.3) 


where  (i  ,j  )  are  the  coordinates  of  the  weight  center, j'  is  the  given  j  coordinate, 
c  c 

and  a,  i  =  0,  ,3  is  the  parameter  of  vector  A,  with  which  we  aproximate  the 

I 

real  ED  curve  Then  we  can  assign 

if  h-'ciiVj . 

IJ  0  otherwise, 

(4.4) 


If  the  approximation  of  ED  curve  in  (3)  is  perfect,  we  can  have 


Q.  -  l  x.  .  ,  i»1 ,  ,m 

l  =i 


and 


l 


Q 


m 


l  iQi-Q,  I. 


(4.5) 


where  Eq  equals  to  zero 

From  the  same  approach  for  P  ,  we  can  also  have 


n 


14  6) 

,and  Ep  equals  to  zero  if  the  estimation  of  ED  curve  is  correct  The  regularity 
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Figure  4-1: 
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constraint  in  our  algorithm  is  much  less  restrictive  than  the  convexity  constraint, 

since  even  the  concave  cross  section  can  also  be  regular.  Therefore,  the  regularity 

constraint  fits  the  ventricular  shapes  much  better.  From  a  normal  shape  of  ventricle, 

which  does  not  have  abrupt  changes  on  its  boundary  to  cause  large  discontinuities 

on  ED  curves,  we  usually  can  expect  Ep  and/or  E^  close  to  zero.  So  the  cross 

section  can  be  reconstructed  if  A  can  be  obtained.  Now  we  state  that  our  problem 

is  to  find  a  optimal  vector  A  which  will  minimize  E  or  E 

P  Q 

In  the  sequel,  let  us  summarize  our  algorithm  of  reconstruction  in  several  steps  as 
follows.  The  detail  of  each  step  will  be  illustrated  in  the  subsequent  subsections 
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Algorithm: 

1.  Find  the  weight  center  (i  ,  j  )  of  two  projection  profiles  by  using 

c  c 

equation  (4.2) 

2.  Assume  ED  is  linear  and  passing  through  (i  ,j  ),  i.e., 

”  C  C 

a 

ED  :  i  -  i  =  a  •  (j-j  ) 

c  1  c 

A 

or  j-j  =  a  •  (i-i  ). 

c  1  c 

By  using  the  coarse-to-fine  approach,  obtain  a^s,  which  minimize  Eq 


3.  Perform  the  consistency  test  to  obtain  the  unique  solution  of  a^  on 
each  layer  of  cross  sections. 

4  Let  ED  be  a  curve  approximated  by  a  third  order  polynomial,  such  that 
ED  :  i-i  =  (a  ±Aa)  +  (a  ±  Aa  Mj-j  )+(a„  ±  Aa  )(j-j  )2+ 

c  0  0  i  i  e  2  2  c 


la3i4a3"ric|3 


(Similarly  for  j  coordinates) 


to  minimize  E  with  respect  to  A=[a  ,a  ,a  ,a  ]'  by  a  gradual  descent 
Q  0  12  3 

method  with  a  fixed  step  AA  The  a*  obtained  in  the  last  step  is  used 
as  the  initial  value  of  a  . 


5.  A  contour  relaxation  procedure  is  then  applied  to  obtain  the  higher 


20 


4.1  A  COARSE-TO-FIIME  APPROACH  TO  FIND  A  * 

i 

For  a  convex  or  slightly  concave  cross  section.  ED  curve  can  roughly  be 
approximated  by  a  line.  The  direction  of  the  estimated  ED  line  implies  the  major 
orientation  of  the  cross  section.  To  obtain  it,  let  us  assume 

i  (  j’  )=  j  +  a  Cj'-j  )  at  j=j ‘ 

C  1  C 

(or  similarly  j  (  i'  )  =  j  +a  (i'-i  ),  at  i=i'  ) 

c  1  c 

(4.7) 

*  * 

and  find  an  optimal  a  such  that  from  the  reconstructed  x  e  {0,1},  we  can  obtain 

ft  • 

E  or  E  .where 
Q  Q 


orE*=Min En=Min 7  IP.-P.I. 

P  a1  Q  a  1  j  J  J 

The  procedure  to  find  E  (Ep)  has  been  described  in  the  previous  subsection  by 
equations  (4.4),(4.5),(4.6). 


Since  the  minimum  of  EQ(or  Ep)  can  not  be  obtained  analytically,  we  use  a  coarse- 

ft 

to- fine  approach  in  finding  a^  As  illustrated  in  Figure  4-2,  we  start  with  an  array 
of  five  equally  spaced  ai  values  covering  the  possible  range  of  ai  [~a,a]  The  three 
neighboring  a^s  which  contribute  the  least  sum  of  three  neighboring  Eqs  (or  Eps) 
are  picked  up  after  each  iteration.  By  inserting  two  new  a^s  into  the  two  slots 
between  three  picked  a^s,  another  array  of  five  equally  spaced  a^s  is  formed  again. 

The  array  is  then  used  to  initiate  another  iteration.  Ultimately,  the  optimal  a*  is 
picked  up  from  one  of  the  five  remaining  a  ’s  which  generate  the  smallest  values  of 

Eq  (or  Ep)  after  four  iterations  of  this  procedure.  Nevertheless,  the  a*  is  not  always 

ft 

unique.  Two  possible  a  ’s  can  be  generated  if  the  project  profiles  of  the  ruler  is 
symmetric  with  respect  to  its  weight  center  The  ED  lines  of  the  reconstructed 

cross  sections  are  symmetrical  with  respect  to  the  axis  passing  through  the  weight 

« 

center  and  parallel  to  the  constructor  In  other  words,  these  two  a  s  are  similar  in 


I  I  is  tli e  selected  center  of  group 
at  each  iteration. 

Figure  4-2: 
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value  but  with  different  signs.  Therefore,  the  other  possible  a  's,  if  available,  can  be 

found  from  the  opposite  portrait  of  the  last  five  remaining  a  's.  In  this 

*  *'  # 
implementation  both  a  and  a  are  recorded  for  the  next  process,  if  the  ai  s  is  not 

larger  than  a  by  15%. 

The  reconstruction  results  from  a  vertica1  and  horizontal  constructor  are  usually 

different  Thus,  the  above  process  can  be  applied  to  both  horizontal  and  vertical 

* 

directions,  but  generate  only  one  set  of  a  ^  from  one  of  the  directions  However, 

both  sets  of  a  are  stored  if  both  directions  generate  similar  E  and  E  This 

1  P  Q 

ambiguity  may  be  caused  by  a  cross  section  highly  symmetric  with  respect  to  the 
origin  The  consistency  test  procedure  in  the  next  subsection  can  be  utilized  to  cope 
with  this  problem. 


4.2.  TEST  OF  CONSISTENCY 

From  a  symmetrical  ruler  projection  profile,  two  cross  sections,  which  are 
symmetrical  to  each  other  with  respect  to  the  direction  parallel  to  the  direction  of 
the  constructor,  can  be  reconstructed  equally  well  Therefore,  the  ambiguity  arises 
when  more  than  one  cross  section  can  be  reconstructed  from  the  same  pair  of 
projection  profiles. 

The  first  step  of  the  reconstruction  process  under  the  assumption  of  a  linear  ED 
is  performed  exclusively  with  the  two  projection  profiles  given  in  the  corresponding 
layer  In  fact,  it  estimates  the  major  orientation  of  the  cross  section  on  the  current 
layer  Thus,  the  major  orientation  of  the  cross  sections  on  all  the  layers  can  be 
obtained  independently  The  orientations  and  error  indices  of  both  possible  answers 
on  every  ambiguous  layer  are  stored  After  computing  cross  sections  on  all  the 
layers  the  ambiguities  can  then  be  resolved  by  using  the  property  that  the  cross 
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sections  on  neighboring  layers  are  smoothly  and  continuously  connected.  Therefore 
their  major  orientations  are  consistent  in  similar  directions. 

To  pick  up  the  real  orientation  from  the  possible  solutions  of  an  ambiguous  layer, 
the  following  two  simple  rules  need  to  be  applied. 

Rulel:  If  both  neighboring  layers  are  unambiguous,  then  pick  the  one  that  gives 
minimum  difference  of  orientations  to  both  neighboring  layers. 

Rule2:  If  only  one  of  the  neighboring  layers  is  unambiguous,  then  pick  the  one 
that  gives  minimum  difference  of  orientation  to  the  unambiguous  cross  section. 

In  general,  the  ambiguities  generated  by  symmetric  projection  profiles  on  the  ruler 
can  be  determined  by  inference  from  the  neighboring  layers. 

The  consistency  of  the  major  orientations  affects  the  3-D  shape  reconstruction 
substantially  only  when  the  reconstructed  cross  sections  are  elongated  in  shape. 
Due  to  the  nature  of  the  algorithm,  the  cross  section  reconstruction  is  more 
accurate  when  ED  lies  along  the  elongated  direction.  Therefore,  the  optimal  ED  lines 
obtained  in  step  one  usually  are  from  the  same  side  of  constructors.  One  from  the 
two  possible  solutions  can  be  determined  easily  by  using  the  two  given  rules.  If  the 
orientation  of  the  cross  section  is  tilted  about  ±45  degrees,  the  reconstructions  by 
using  horizontal  or  vertical  projection  profiles  as  constructors  can  both  obtain  fairly 
good  results.  It  is  preferred  to  use  the  same  side  of  constructor  as  its  two 
neighboring  layers.  However,  the  two  neighboring  layers  may  not  use  the  same  side 
of  constructors  in  the  transition  area  where  the  orientations  of  the  cross  sections 
change  from  less  than  to  more  than  ±45  degrees.  The  two  given  rules  are  still 
applicable  since  the  real  ED  lines  of  the  horizontal  and  the  vertical  constructors  are 
in  similar  orientations  when  their  major  direction  of  the  cross  section  are  around 
±45  degrees.  The  orientation  that  generates  minimum  projection  error(Ep  or  Eq)  and 
satisfies  the  two  given  rules  is  picked  up  in  the  transition  case. 

The  process  of  the  consistency  test  can  be  performed  from  top  layer  to  bottom 
layer  (and  then  in  reverse  if  necessary!  by  using  the  results  from  the  first  step 
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It  is  comparatively  less  expensive  in  computation  because  point  to  point  matching 
is  unnecessary  in  our  approach.  Instead,  the  properties  of  ED  line  can  be  utilized 
Thus  ,  the  consistent  shape  of  3-D  objects  is  available,  after  this  test,  for  the 
detailed  adjustment  in  the  following  steps. 


4.3  ITERATIVE  DESCENT  PROCEDURE  TO  ESTIMATE  ED 
CURVE  BY  A  POLYNOMIAL 

The  ED  curve  is  not  linear  on  some  of  the  convex  cross  sections  and  most  of 
the  concave  cross  sections.  The  ED  line  obtained  from  the  above  process  needs  to 
be  adjusted  to  better  fit  the  real  ED  curve.  Therefore,  a  third  order  polynomial  is 
used  to  approximate  the  ED  curve  under  the  constraint  that  the  weight  center  of 
the  cross  section  is  invariant  A  higher  order  polynomial  is  intuitively  easier  to  fit  an 
arbitrary  curve  In  practice,  a  third  order  polynomial  is  usually  moderate  to  simulate  a 
nonlinear  ED  curve  and  computationally  less  expensive 

It  is  known  that  the  weight  center  of  the  real  cross  section  can  be  uniquely 
determined  from  the  given  two  projection  profiles.  So.  the  weight  center  of  the 
reconstructed  cross  section  certainly  should  be  always  kept  identical  to  the  real  one 
As  a  result  of  this  constraint,  only  three  parameters  need  to  be  adjusted  in  the 
optimization  process.  Now,  we  begin  to  further  refine  the  ED  curve  by  using 

i  ( j  ' )  ■  i  c+  (aQ®+Aag)+  (a^tAa^)  (j  ' -jc)  +  (a^± Aaj)  (j  ' -j c)  2+ 

(a^  +  Aa^)  (j  '  -jc)  3 

148) 

(Similarly,  for  the  other  projection 

j  (i  ' )  ■  j c+  (aQ®+AaQ)+  (a^  +  Aaj)  (i'-'c)  +  (a2°±Aa2)  ( i  '  -  i  c)  2+ 

(aj°±Aaj)  (i  1  -ic)  3 

) 
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At  first,  let  a  be  the  initial  value  of  a  and  zero  be  the  initial  values  of  a  _  a  ,  a 
1  1  0  2  3 

to  form  an  initial  vector  of  parameters,  denoted  by  Aq  A  fixed  array  of  Aa , 
i=  1,2,3  are  given  such  that  the  external  boundary  points  are  changed  by  at  least  one 

cell  when  the  corresponding  a  varies  with  the  amount  Aa.  The  change  of  a  value 

i  i  o 

simply  moves  the  weight-  center  of  the  reconstructed  cross  section  away  from  the 
real  weight  center  So  the  aQ  value  is  only  adjusted  by  the  other  parameters  in 
order  to  compensate  the  generated  deviations  to  the  weight  center  when  those 
parameters  are  changed.  The  iterative  descent  procedure  is  therefore  applied  only 
to  three  parameters  of  A^;  aQ  is  kept  updated  when  necessary.  Four  points 
1  0  0 

denoted  as  a  =  a  ±  Aa ,  where  a  has  not  been  adjusted  in  the  previous  iteration 

i  i  i  i 

# 

and  i  e  {1,2,3},  are  checked  to  find  the  optimal  a  which  results  in  a  minimum 

I 

E  (E  ),  called  E  (E  ).  If  E  (E  )  <  E°(E°),  then  a°  =  a*  This  process  is  repeated  till  no 
QP  QP  QPQP  li 

further  improvement  can  be  made  The  final  vector  of  A°  is  then  the  optimal 
solution  of  the  problem  that  generates  minimum  sum  of  error  between  the  real  and 
estimated  projection  profiles 


4.4  CONTOUR  RELAXATION 

The  real  ED  curve  of  a  cross  section  can  fairly  be  approximated  by  a  third  order 
polynomial.  However,  the  attempt  to  obtain  the  exact  ED  curve  still  involves  some 
higher  order  terms.  To  resolve  these  higher  order  terms  by  increasing  the  order  of 
the  polynomial  which  is  used  to  approximate  ED  curve  is  too  computationally 
expensive  The  contour  relaxation  algorithm  is  designed  to  solve  this  problem 
without  directly  computing  the  higher  order  polynomials  Based  on  the  theory,  it 
minimizes  the  difference  between  the  given  and  estimated  projection  profiles  in 
order  to  minimize  the  difference  between  given  and  estimated  cross  sections  under 
the  condition  that  the  regularity  constraint  on  the  contour  of  the  cross  section 
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should  not  be  violated.  The  difference  between  the  real  and  estimated  ED  curves 
introduces  error  represented  as  the  difference  between  the  estimated  and  given 
|  projections  along  the  ruler  coordinate.  It  is  clear  from  our  proposition  that  the  real 

;  cross  section  can  be  uniquely  reconstructed  from  a  given  constructor  and  its 

i  corresponding  real  ED  curve.  Therefore,  the  contour  relaxation  algorithm  intends  to 

|  eliminate  the  difference  between  given  and  estimated  projection  profiles  to  obtain 

the  estimated  ED  as  correct  as  possible.  Now  let  us  introduce  the  algorithm  in  the 
following. 

|  To  simplify  the  illustration,  we  use  a  horizontal  constructor  as  in  Figure  4-3;  the 

'  alternative  case  can  be  solved  similarly.  First,  let  us  set  the  projection  error  e(j)  as 

I  the  projection  of  estimated  cross  section  minus  the  real  given  projection  at  each 

row  j  along  the  ruler  coordinate.  The  upward  potential  U(k),  and  downward  potential 

i  D(k),  which  will  be  explicitly  defined  later,  at  each  column  i  along  the  constructor 

j  i 

i  coordinate  are  used  as  the  basic  measures  in  this  algorithm.  As  we  know,  the 

!  deviation  between  real  and  estimated  ED's  generates  the  difference  between 

estimated  and  ruler  projection  profiles,  and  contributes  error  to  the  projection  error 
e(j).  This  error  can  obviously  be  compensated  by  shifting  the  estimated  ED  upward 
or  downward  toward  the  real  ED.  However,  the  real  ED  is  unknown.  Therefore,  the 
adopted  U(k)  and  D  (k)  represent  the  potential  to  reduce  error  by  shifting  upward  or 
downward  along  column  i  from  point  k- 1  to  point  k  away  from  the  current 
estimated  ED.  It  is  obvious  that  the  potential  of  shifting  upward  or  downward 
increases  when  large  projection  error  can  be  compensated  after  the  shifting.  From 
Figure  4-3,  we  define  U  (k)  and  D(k)  as  below 

i  i 


Let  m  and  n  denotes  the  row  coordinates  of  the  upper  and  lower  boundaries  of 

i  i 

the  estimated  cross  section  at  column  i.  Moving  this  column  one  point  upward  (k=  1 ) 
will  make  the  projection  error  e(m  +  1)  increase  one  while  e(n)  will  decrease  one. 

I  I 

Therefore,  larger  e(n)  and  smaller  elm  +  1)  implies  more  potential  to  move  this 

I  I 

column  upward.  Thus,  we  define 
U  ( 1 )  =  e(n  )  -  elm  + 1 ), 

i  i  i 


and 
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U  (k)  =  e(n+k-1)  -  e(m +k), 

I  i  I 

when  the  column  is  moved  from  point  k- 1  to  point  k  upward.  Similarly,  for 
shifting  downward,  we  define 

D(1)  =  e(m)  -  e(n  —  1 ), 

i  i  i 

and 

D  (k)  =  e(m-k+1)  -  e(n -k). 

i  i  i 

Using  U  and  D,  we  then  define 

TU  =  2  •  U(1)  +  U<2)  +  U  (3) 

i  i  i  i 

TD  =2  •  D(1)  +  D  (2)  +  D  (3) 

i  jii 

and 

f(i)  =  TU,  d(i)=1,  if  TU  -  TD  >  T 

i  i  i  i 

f(i)  =  TD,  d(i)=- 1,  if  TD  -  TU  >  T, 

i  iii 

f(i)  =  0,  d(i)=0,  ortherwise, 

where  T  =  0.75  Max(  TU  ,  TD  ). 

i  i  i 

To  include  the  influence  from  the  neighboring  columns,  we  then  assign 
g(i>  =  d(i-1)*f(i-1)  +  2d(i)  •  f  (i>  +  d(i+1)*f(i+1), 

to  be  used  in  the  iterative  process. 

The  TU  (or  TD)  from  the  linear  combination  of  U  (or  D)  assigns  more  weight  to 

II  II 

the  columns  need  to  be  shifted  more  deeply  upward(downward).  The  f (i)  and  d(i)  are 
then  defined  so  that  each  column  i  can  only  be  shifted  upward,  downward,  or  kept 
still.  The  larger  the  f(i)  implies  the  more  potential  to  move  this  column.  Considering 
neighboring  columns  are  usually  shifted  smoothly  in  the  same  direction,  g(i)  grants 
priority  to  the  column  where  both  its  neighbors  and  itself  have  a  strong  tendency  to 
be  adjusted.  Now  let  us  state  the  contour  relaxation  algorithm  as  follows: 


1.  Adjust  those  columns  which  exceed  the  upper  and  lower  end  points  of 
the  ruler  to  assure  the  contour  is  within  the  possible  region  of  cross 
section 

2.  Let  k=0  Obtain  the  initial  e  (j)  for  each  row  on  the  ruler  from  the 
estimated  cross  section. 

k.  k  k- 1 

3.  k=k+ 1 .  Compute  g  (i)  and  d  (i)  from  e  (i)  for  each  column  i  on  the 
constructor 

k-  1 

4.  For  each  row  j  with  e  (j)  P  0,  pick  up  a  column  with  maximum 

k  k  *  ^ 

g  (i)  among  those  candidates,  which  reduce  |  e  (j)  |  by  shifting 

k 

according  to  d  (i),  without  violating  the  regularity  constraint  of  the 
boundary 

5.  For  each  picked  column  i,  adjust  it  by  shifting  toward  the  direction 

k  k  k 

given  by  d  (j)  and  then  set  d  (j)=0  If  the  d  (j)  has  been  set  to  zero,  do 

k 

no  adjustment  Compute  e  (i)  after  all  columns  have  been  adjusted 


|r 

6.  If  E=  I  |  e  (j)  |  is  reduced,  go  to  2,  else  stop 


Figure  4-4  demonstrates  an  example  for  the  use  of  the  contour  relaxation 
algorithm.  Figure  4-4(a)  is  obtained  from  the  ED  curve  approximated  by  a  third 
order  polynomial.  Figure  4— 4(b)  shows  the  results  after  applying  this  algorithm 
Considerable  improvement  has  been  made  by  the  contour  relaxation  In  fact,  the 
algorithm  can  generally  further  refine  the  reconstructed  cross  section  as  long  as  the 
regularity  is  satisfied  The  higher  order  terms  of  the  estimated  ED  curve  are 
resolved  implicitly 


4.5  SOME  EXAMPLES  OF  THE  CROSS  SECTION 
RECONSTRUCTION 


Some  examples  of  the  cross  section  reconstruction  are  given  in  this  section 
Figures  4-5  to  4-10  illustrate  some  of  the  experimental  results  with  our  algorithm 
for  six  randomly  drawn  cross  sections  Each  cross  section  is  drawn  by  a  graphic 
program  into  a  60  x  60  matrix  The  projection  profiles  are  constructed  by  summing 
up  the  elements  belonging  to  the  region  of  cross  section  columnwise  or  rowwise 
The  reconstructed  cross  section  is  compared  with  the  original  The  relative  mean 
error  used  in  (6,1  1,12)  is  adopted  to  measure  the  performance  It  is  defined  as 


The  Q  are  the  elements  of  the  original  cross  section  and  x  are  the  elements  of 
ij  ij 

the  reconstructed  one  As  mentioned  above,  the  reconstruction  does  not  refer  to 
any  predefined  binary  mask  models  However  the  results  are  quite  good  even  when 
the  original  cross  section  is  very  irregular  Table  4-1  summarizes  the  experimental 
results  for  the  cross  sections  shown  in  Figure  4-5  to  4-10.  For  each  figures  set. 
(a)  shows  the  original  cross  section,  (b)  illustrates  the  reconstructed  results  and  the 
estimated  ED  curve  given  by  our  method,  (c)  is  the  reconstruction  results  obtained 
by  the  ellipse  approximation  method.  The  dark  curve  inside  the  reconstructed  region 
in  (b)  is  the  ED  curve  obtained  by  our  algorithm.  From  50  randomly  drawn  pictures, 
the  reconstruction  results  for  a  regular  shape  ususlly  can  achieve  more  than  96% 
conformity,  which  is  better  than  th*  results  shown  in  (4  5).  The  conformity  measure 
is  given  by  the  number  of  elements  common  to  the  original  and  the  reconstructed 
cross  sections  Thus  the  conformity  measure  is  one  hundred  minus  half  of  the 
relative  error  R  in  percent  As  to  the  cross  sections  of  irregular  shape,  we  usually 
get  less  than  20%  R.  in  other  words,  it  is  around  90%  conformity  Figure  4-5  to 
4-10  show  fairly  regular  cross  sections:  the  results  are  as  good  as  predicted 
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Figure  4-9  and  4- 1 0  show  cross  section  with  obvious  irregularity,  which  is 
generated  by  large  concavity  and  sharp  changes  on  boundary.  Although  the  cross 
sections  are  irregular,  the  reconstructed  results  still  keep  the  approximate  shapes 
for  visual  recognition.  In  the  next  section,  we  are  going  to  do  the  experiments  with 
real  X-ray  pictures 


Table  4-1:  Experimental  results  for  Figures  4-5  to  4-10 


Figure  Total  No.  Our  Method  Ellipse  Approximation 

No.  of  Points 


CHAPTER  5 


AN  EXPERIMENTAL  RESULT 

In  this  section,  we  will  demonstrate  an  example  with  two  X-ray  films  of  a  bag  of 
dye  to  show  the  applicability  of  the  object  reconstruction  technique  described  in 
Section  3.  The  procedure  of  object  reconstruction  is  shown  in  Figure  5-1.  The 
films  were  digitized  by  an  Optronix  drum  scanner  into  512  by  512  arrays  of  pixels. 
The  grey  value  of  each  pixel  is  recalled  by  an  8-bit  integer  which  scales  the 
intensity  of  brightness  into  the  range  between  0  and  255.  Low  grey  value 
represents  low  intensity  (dark  pixels),  while  high  grey  value  represents  high  intensity 
(bright  pixels).  The  second  step  is  to  apply  the  logarithmic  transformation  to  both 
left  and  right  digitized  images,  such  that  the  exponential  absorption  of  the  radiation 
is  compensated.  Then,  the  grey  value  of  the  dye  free  region  is  subtracted  from  the 
image  since  we  are  only  interested  in  the  net  X-ray  absorption  by  dye,  which  is 
proportional  to  the  depth  of  the  ventricle. 

For  a  real  ventricular  reconstruction,  we  have  to  do  some  averaging  work  to 
assure  the  quality  of  image  because  the  mixture  of  the  injected  dye  is  usually 
incomplete.  However,  this  experiment  utilizes  a  full  bag  of  dye,  so  the  incomplete 
mixture  problem  does  not  have  to  be  considered  in  this  case.  The  next  step  is  the 
ventricular  boundary  detection  Several  papers  have  been  presented  for  this  problem. 
Pope  et  al  utilize  the  Dynamic  Search  Algorithm  (11)  and  Bocker  utilizes  the 
Laplacian-Gaussian  operator  ( 1 2),  which  are  pretty  reliable  for  the  ventricular 
boundary  detection  .  In  this  experiment,  we  utilize  a  simple  threshold  method  to 
segment  the  region  of  object  The  detected  boundary  is  satisfactory  since  the  used 
object  is  150  cc  of  pure  dye  in  a  plastic  bag  whose  shape  is  regular  and 
unambiguous 


digitization 
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After  the  boundary  detection,  the  segmented  region  is  smoothed  by  a  medium 
filter  for  eliminating  the  corrupted  noise  at  the  present  experiment  The  processed 
images  are  ready  to  use  except  their  sizes  are  different  since  they  were  taken  by 
two  different  sets  of  cameras.  Therefore,  these  two  images  are  scaled  and  aligned 
so  that  we  can  slice  the  two  images  into  corresponding  arrays  of  pixels  Each  cross 
section  will  then  be  reconstructed  from  the  grey  values  of  the  corresponding  pair 
of  pixel  arrays.  However,  we  can  not  apply  these  grey  values  directly  into  the 
binary  reconstruction  process  because  grey  value  is  not  equivalent  to  the  depth  of 
the  object  Fortunately,  the  depth  of  the  object  and  grey  level  are  proportional  to 
each  other.  It  can  be  expressed  by 

D  (i)  =  a* I  (i) 

X  X 

where  D  is  the  depth,  I  is  the  grey  value,  and  a  is  an  adjustment  coefficient 

X  X 

decided  by  the  characteristics  of  the  X-ray  machine.  In  this  experiment,  the  a 
values  for  both  projections  are  unknown.  For  simplicity,  we  assume  the  highest  grey 
value  corresponds  to  the  width  of  the  grey  level  profile  on  the  other  projection  So 
these  two  adjustment  coefficients  can  be  decided  and  used  to  transfer  grey  level 
profiles  into  depth  profiles.  The  generated  depth  profiles  are  then  used  to 
reconstruct  the  cross  sections.  The  reconstructed  cross  sections  are  put  into  a 
stack  to  represent  the  three-dimensional  objects.  On  the  other  hand,  the  volume  of 
the  object  is  also  computed  with  a  Simpson  integration  formula. 

Volume  =  h/3  (  Z  4  •  area)  +  (  Z  2  •  area) 

odd  «ven 

Figure  5-2  shows  the  original  X-ray  pictures  for  both  left  and  right  projections. 
Figure  5-3  shows  the  same  image  after  preprocessing.  Figure  5-4  shows  the 
reconstructed  three-dimensional  object  by  piling  up  each  reconstructed  cross 
section.  The  estimated  volume  of  this  object  is  153.2  cc,  which  is  around  97.9% 
accurate  with  respect  to  the  real  volume  150  cc.  The  ellipse  method  by  the  HP. 
ventricle  estimation  system  obtains  1 63  cc  for  the  object,  which  is  around  9 1 .4% 
accurate.  Our  method  obtains  better  volume  estimation,  and  the  reconstructed  shape 
(Fig.  5-4)  is  close  to  the  real  shape  of  the  object  according  to  doctor's  observation. 


CHAPTER  6 


CONCLUSIONS  AND  SUGGESTIONS 
FOR  FURTHER  IMPROVEMENT 

In  this  report,  we  present  a  new  method  to  reconstruct  the  three-dimensional 
shape  and  estimate  the  volume  of  an  object  by  using  its  two  orthogonal  projections. 
We  have  proved  that  the  cross  section  can  be  optimally  reconstructed  by  minimizing 
an  error  index  which  is  generated  by  the  sum  of  the  absolute  difference  between 
the  original  and  the  estimated  projection  profiles.  This  error  index  is  uniquely  and 
continuously  proportional  to  the  sum  of  error  between  the  original  and  the  estimated 
cross  sections  under  the  condition  that  the  original  cross  section  is  regular  and  has 
a  monotonically  non-decreasing  or  non-increasing  equal-divisor  curves.  A  new 
algorithm  has  been  designed  based  on  these  properties. 

The  algorithm  has  been  tested  on  synthetic  cross  section  as  well  as  real  X-ray 
pictures.  From  synthetic  pictures,  it  achieves  better  than  96%  conformity  for  regular 
cross  sections  and  better  than  91%  conformity  for  irregular  ones.  The 
reconstruction  for  X-ray  pictures  is  also  better  than  95%  conformity  in  every  layer 
of  cross  section.  The  processing  speed  is  less  than  3  seconds  for  each  cross 
section  by  a  Fortran  implementation  in  VAX  780.  Thus,  the  algorithm  are  reliable  and 
fast  for  the  different  experiments  that  we  have  done  so  far. 

In  the  next  phase  of  this  project,  we  are  going  to  apply  this  algorithm  to  real 
ventricular  images  To  obtain  better  results  in  real  ventricular  images,  we  suggest  to 
further  investigate  the  adjustment  coefficient  a  of  X-ray  absorption.  We  also 
propose  some  methods  to  depress  the  corrupted  noise  in  real  X-ray  images.  At 
last,  some  more  detail  consideration  will  be  made  to  directly  reconstruct  3-D 
objects  from  the  two  projection  images  Some  more  studies  for  the  irregular  case 
is  also  of  interest 
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